library(adehabitat)
library(fields)
library(maptools)

#set the working directory
workdir = "C:/working/maxent_RFspp_3to6/Maxent/output_14Mar2008/"
setwd(workdir)

#get the subregion mask
subregion = read.shape('C:/working/maxent_RFspp_3to6/Maxent/enviro_data/subregions_new07.shp')
subregion = Map2poly (subregion)

#get the occurrences
occur = read.csv("C:/working/maxent_RFspp_3to6/Occur_14Mar2008.csv")
occur = data.frame(spp=occur$species_ID,x=occur$east,y=occur$north)

#get the files in the potential and realized distributions
pot.spp = list.files(paste(workdir,"potential_dist/",sep=""))
t.spp = c(NA) #need to remove any files from list that are not ascii grid files
for (pspp in pot.spp) {
  if (length(grep(".asc",pspp))>0) t.spp = c(t.spp,pspp)
}
pot.spp = na.omit(t.spp); rm(t.spp)
real.spp = list.files(paste(workdir,"realized_dist/",sep=""))

#check if there is both a realized and potential distribution
for (pspp in pot.spp) {
  spp = gsub(".asc","",pspp) #get the spp id
  xy = data.frame(x=occur$x[occur$spp==spp], y=occur$y[occur$spp==spp]) #get the occurrences for the species
  print(spp)
  both = F
  for (rspp in real.spp) {
    if (rspp == pspp) both = T
  }
  if (both==T) { #if both realized and potential distributions exist
    png(paste("dist_plots/",spp,".png",sep=""),width=1000,height=1000)
      par(mfrow=c(1,2),pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      tasc = import.asc(paste("potential_dist/",spp,".asc",sep=""))
      image(tasc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c("grey",heat.colors(5)[seq(5,1)])) #plot the raster
      plot.polylist(subregion,add=T,border="grey30") #add the subregions
      points(xy,pch='+')  #add the occurrences

      tasc = import.asc(paste("realized_dist/",spp,".asc",sep=""))
      image(tasc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c("grey",heat.colors(5)[seq(5,1)]))  #plot the raster
      plot.polylist(subregion,add=T,border="grey30") #add the subregions

      mtext("            Potential Distribution                         Realized Distribution", line = 1.5, outer = TRUE, side=1, cex = 2.1)
      mtext(spp, line = 1.5, outer = TRUE, side=3, cex = 2)
    dev.off()
  } else { #if only the potential distribution exists
    png(paste("dist_plots/",spp,".png",sep=""),width=500,height=1000)
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      tasc = import.asc(paste("potential_dist/",spp,".asc",sep=""))
      image(tasc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c("grey",heat.colors(5)[seq(5,1)]))  #plot the raster
      plot.polylist(subregion,add=T,border="grey30") #add the subregions
      points(xy,pch='+')  #add the occurrences
      
      mtext("Potential Distribution", line = 1.5, outer = TRUE, side=1, cex = 2.1)
      mtext(spp, line = 1.5, outer = TRUE, side=3, cex = 2)
    dev.off()
  }
}

